Рассмотрим временной ряд динамики средней цены (доллар/фунт) на грейпфруты в США в Цетрально-Западном и Северо-Восточном регионе за двадцатитрехтилетний период (длина ряда 276 точек) и его периодограмму.
data.mw <- read.csv('APU0200711411.csv', sep = ',', header = T)
data.ne <- read.csv('APU0100711411.csv', sep = ',', header = T)
df <- data.frame(data.mw$APU0200711411[1:276], data.ne$APU0100711411[1:276])
colnames(df) <- c("Midwest", "Northeast")
price <- ts(df, start = c(1980, 1), end = c(2002,12), frequency = 12)
plot(price)
spec.pgram(price, log='no', detrend = FALSE, taper = 0)
Графики и периодограммы очень схожи. Неудивительно, ведь это цена на одни и те же фрукты в разных частях страны. Так как сигнал условно одинаковый, дополнительная информация из второго канала (регион Northeast) может помочь лучше отделиться от шума. Сравним разложение цены Midwest методом SSA и M-SSA.
s.mw <- ssa(price[,1])
plot(wcor(s.mw, groups = 1:20))
plot(s.mw, type = "vectors", idx = 1:15)
plot(s.mw, type = "paired", idx = 1:15)
Видим, что 2-3 компонента перемешалась с трендом.
Попытаемся улучшить разделимость внутри сигнальных компонент с помощью DerivSSA.
fs.mw <- fossa(s.mw, nested.groups = c(1:7))
plot(wcor(fs.mw, groups = 1:20))
plot(fs.mw, type = "vectors", idx = 1:15)
plot(fs.mw, type = "paired", idx = 1:15)
Разделимость стала намного лучше.
Посмотрим на восстановленный сигнал и acf остатков.
r.mw <- reconstruct(fs.mw, groups = list(Trend = c(3:7),
Seasonality = c(1:2)))
plot(r.mw, add.residuals = FALSE)
acf(residuals(r.mw))
По acf остатков можем предположить, что это шум.
Теперь применим Multi-channel SSA.
s <- ssa(price, kind = "mssa")
plot(wcor(s, groups = 1:20))
plot(s, type = "vectors", idx = 1:15)
plot(s, type = "paired", idx = 1:15)
Компоненты, в том числе полугодовая гармоника, хорошо отделились от шума, тренд не перемешан с гармониками сигнала.
Посмотрим на восстановление и сравним полученные разложения.
r <- reconstruct(s, groups = list(Trend = c(1, 4, 5, 8, 9),
Seasonality = c(2:3, 6:7)))
plot(r, add.residuals = FALSE,
plot.method = "xyplot",
superpose = TRUE, auto.key = list(columns = 3))
plot(r$Trend[,1], col = "green", ylab = "Trend (Midwest)")
lines(r.mw$Trend, col = "red")
legend("topleft",
legend = c("MSSA", "DerivSSA"),
col = c("green", "red"), lty = 1)
plot(r$Seasonality[,1], col = "green", ylab = "Seasonality (Midwest)", ylim = c(-0.15, 0.17))
lines(r.mw$Seasonality, col = "red")
legend("topleft",
legend = c("MSSA", "DerivSSA"),
col = c("green", "red"), lty = 1)
spec.mw <- spec.pgram(residuals(r.mw), log='no', detrend = FALSE, plot = FALSE)
spec <- spec.pgram(residuals(r)[,1], log='no', detrend = FALSE, plot = FALSE)
plot(residuals(r.mw), col = "red", ylab = "Residuals (Midwest)", ylim = c(-0.15, 0.17))
lines(residuals(r)[,1], col = "green")
legend("topleft",
legend = c("MSSA", "DerivSSA"),
col = c("green", "red"), lty = 1)
plot(spec.mw$freq, spec.mw$spec, col = "red", type = "l",
xlab = "frequency", ylab = "spectrum", main = "Периодограмма остатков")
lines(spec$freq, spec$spec, col = "green",)
legend("topright",
legend = c("MSSA", "DerivSSA"),
col = c("green", "red"), lty = 1)
Результаты схожи, однако применение M-SSA позволило нам отделить компоненты полугодовой сезонности от шума.
Загрузим изображение, добавим к нему текстуру и посмотрим, сможет ли 2D-SSA выделить изображение.
raw <- readRaw('image.raw')
matrix <- matrix(as.numeric(raw$fileRaw), nrow = 400, ncol = 400, byrow = TRUE)
x <- 1:(400*400)
noise <- matrix(5*sin(x)+10*cos(2*x), nrow = 400, ncol = 400, byrow = TRUE)
image <- matrix + noise
s <- ssa(image, kind = "2d-ssa", L = c(25, 25))
plot(s, type = "vectors", idx = 1:20,
cuts = 255, layout = c(10, 2),
plot.contrib = FALSE)
plot(wcor(s, groups = 1:30), scales = list(at = c(10, 20, 30)))
r <- reconstruct(s, groups = list(Noise = c(7:8, 14:15)))
plot(r, cuts = 255, layout = c(3, 1))
Искусственно наложенный шум хорошо отделился.